#include "SphData.hpp"

SphData::SphData()
{   }

void SphData::clear()
{
    for(int i = 0; i < dnames.size(); ++i){
        dvalues[i].clear();
    }
    dnames.clear();
    map_names.clear();
}
    
void SphData::readFile(std::string s)
{
    int i, j, n;
    double t0, t1, t2, t3;        
    
    std::string line, sub;
    std::vector<std::string> file;
    std::ifstream infile (s.c_str(), std::ios_base::in);
    
    t0 = get_time();
    
    while (getline(infile, line, '\n'))
    {
        file.push_back(line);
    }
    
    infile.close();

    std::cout << "Read " << file.size() << " lines.\n";
    std::cout << "Time to get file: "<< get_time() - t0 << "\n";
    
    /* time */

    line = file[28]; 
    j = line.find("= ");
    sub = line.substr(j + 2, 12);
    tdata = strtod(sub.c_str(), NULL);
    
    std::cout << "Time = " << tdata << "\n";
    
    /* bounds */
    
    {
        line = file[33];
        j = line.find("= ");
        sub = line.substr(j + 2, 14);
        xLeft = strtod(sub.c_str(), NULL);
        
        line = file[34];
        j = line.find("= ");
        sub = line.substr(j + 2, 14);
        xRight = strtod(sub.c_str(), NULL);
        
        std::cout<<"Left-right: "<<xLeft<<" "<<xRight<<"\n";
        
        line = file[35];
        j = line.find("= ");
        sub = line.substr(j + 2, 14);
        yDown = strtod(sub.c_str(), NULL);
        
        line = file[36];
        j = line.find("= ");
        sub = line.substr(j + 2, 14);
        yUp = strtod(sub.c_str(), NULL);
        
        std::cout<<"Down-up: "<<yDown<<" "<<yUp<<"\n";
    }
        
    /* names */
    
    line = file[38];
    std::istringstream iss(line);
    while(iss){
        iss >> sub;
        if(sub.size()) dnames.push_back(sub);
    }
    
    /* data */
    double value;
    
    t1 = get_time();
    
    dvalues.resize(dnames.size());
    for(j = 0; j < dnames.size(); ++j){
        map_names[dnames[j]] = j;
        dvalues[j].empty();
    }

    std::cout<<"Names:"<<dnames.size()<<"\n";

    double sum = 0.0;
    int jj;
    
    std::set<int> set_names;
    std::vector<std::string> names;
    
    names.push_back("X");
    names.push_back("Y");
    names.push_back("P");
    names.push_back("D");
    names.push_back("UX");
    names.push_back("BURN");

    for(jj = 0; jj < names.size(); ++jj) set_names.insert(map_names[names[jj]]);
    
    int min_j = *min_element(set_names.begin(), set_names.end());
    int max_j = *max_element(set_names.begin(), set_names.end());
    
    t2 = get_time();
    std::cout<<"Time to read names:"<<t2 - t1<<" "<<min_j<<" "<<max_j<<"\n";
    
    int k, k0, l;
    std::string s1;
    for(i = 39; i < file.size(); i++){
        if(dkeys.size() == 0 || dkeys[i - 39]){
            //std::istringstream iss(file[i]);
            std::string s = file[i];
            int l = s.size();
            
            k = 0;
            while(k < l && s[k] == ' ') k++;
            k0 = k;
            
            n = dnames.size();
            for(j = 0; j < n && j <= max_j; j++){
                
                while(k < l && s[k] != ' ') k++;
                
                if(set_names.find(j) != set_names.end()){
                    s1 = s.substr(k0, k-k0);
                    value = strtod(s1.c_str(), NULL);
                    dvalues[j].push_back(value);
                }else{
                    value = 0;
                }
                
                sum = sum + value;
                
                while(k < l && s[k] == ' ') k++;
                k0 = k;
            }
        }
    }
    
    std::cout<<"Sum:"<<sum<<"\n";
    
    t3 = get_time();
    std::cout<<"Time to read:"<<t3 - t2<<"\n";
}

/* compare */

void SphData::readFileCompare(std::string s)
{
    int i, j, n;
    double t0, t1, t2, t3;        
    
    std::string line, sub;
    std::vector<std::string> file;
    std::ifstream infile (s.c_str(), std::ios_base::in);
    
    t0 = get_time();
    
    while (getline(infile, line, '\n'))
    {
        file.push_back(line);
    }
    
    infile.close();

    std::cout << "Read compare " << file.size() << " lines.\n";
    std::cout << "Time to get file compare: "<< get_time() - t0 << "\n";
    
    /* data */
    double value, value1, err_value;
    
    t1 = get_time();

    std::cout<<"Names:"<<dnames.size()<<"\n";

    double sum = 0.0;
    int jj;
    
    std::set<int> set_names;
    std::vector<std::string> names;
    
    names.push_back("X");
    names.push_back("Y");
    names.push_back("P");
    names.push_back("D");
    names.push_back("UX");
    names.push_back("BURN");

    for(jj = 0; jj < names.size(); ++jj) set_names.insert(map_names[names[jj]]);
    
    int min_j = *min_element(set_names.begin(), set_names.end());
    int max_j = *max_element(set_names.begin(), set_names.end());
    
    t2 = get_time();
    std::cout<<"Time to read names:"<<t2 - t1<<" "<<min_j<<" "<<max_j<<"\n";
    
    int k, k0, l;
    int j_value;
    std::string s1;
    for(i = 39; i < file.size(); i++){
        if(dkeys.size() == 0 || dkeys[i - 39]){
            //std::istringstream iss(file[i]);
            std::string s = file[i];
            int l = s.size();
            
            k = 0;
            while(k < l && s[k] == ' ') k++;
            k0 = k;
            
            n = dnames.size();
            j_value = 0;
            for(j = 0; j < n && j <= max_j; j++){
                
                while(k < l && s[k] != ' ') k++;
                
                if(set_names.find(j) != set_names.end()){
                    s1 = s.substr(k0, k-k0);
                    value = strtod(s1.c_str(), NULL);
                    value1 = dvalues[j][j_value];
                    
                    err_value = (abs(value) + abs(value1) > 1.0E-32)?2.0 * abs(value - value1) / (abs(value) + abs(value1)):0;
                    
                    dvalues[j][j_value] = err_value;
                }else{
                    value = 0;
                }
                
                sum = sum + value;
                
                while(k < l && s[k] == ' ') k++;
                k0 = k;
            }
        }
    }
    
    std::cout<<"Sum:"<<sum<<"\n";
    
    t3 = get_time();
    std::cout<<"Time to read:"<<t3 - t2<<"\n";
}

/* data */

std::vector<double> SphData::getData(std::string name)
{
    int j = map_names[name];
    std::vector<double> aa = dvalues[j];
    std::cout<<"Size of <"<<name<<"> is "<<aa.size()<<"\n";
    return aa;
}

void SphData::smoothData(std::vector<double> & xx, std::vector<double> & vv)
{
    std::vector<double> xx2, vv2, ss2;
    
    int i, j;
    int nn = xx.size();
    int nn2 = int (nn / 33.3);
    
    double min_x, max_x;
    min_x = *min_element(xx.begin(), xx.end());
    max_x = *max_element(xx.begin(), xx.end());
    
    xx2.resize(nn2 + 1, 0.0);
    vv2.resize(nn2 + 1, 0.0);
    ss2.resize(nn2 + 1, 0.0);
    
    
    for(j = 0; j < nn2; ++j){
        xx2[j] = min_x + j * (max_x - min_x) / nn2;
    }
    
    for(i = 0; i < nn; ++i){
        j = int(nn2 * (xx[i] - min_x) / (max_x - min_x));
        vv2[j] += vv[i];
        ss2[j] += 1.0;
    }
    
    for(j = 0; j < nn2; ++j){
        if(ss2[j] > 0) vv2[j] = vv2[j] / ss2[j];
    }
    
    xx = xx2;
    vv = vv2;
}

double SphData::detectFront(std::vector<double> & xx, std::vector<double> & yy)
{
    int i, j, n, nn;
    std::vector<std::pair<double, double> > xxyy;
    std::pair<double, double> xy;
    
    nn = xx.size();

    double average_y, sum_y;

    double min_y, max_y;
    min_y = *min_element(yy.begin(), yy.end());
    max_y = *max_element(yy.begin(), yy.end());
    
    /* average */
    
    sum_y = 0.0;            
    for(i = 0; i < nn; ++i) sum_y += yy[i];
    average_y = sum_y / nn;
    
    /* front */
    
    double x0;
    double delta = 0;
    double max_delta = 0;
    
    int i0;
    int ik;
    int ik1 = 1;
    int ik2 = 15;
    
    double y_left, y_right;
    double x_left, x_right;
    
    for(i = ik1 + ik2; i < xx.size() - ik1 - ik2; ++i){
        delta = 0;
        
        /* up */
        
        y_left = 0;
        x_left = 0;
        for(ik = i - ik1 - ik2; ik < i - ik1; ++ik){
            y_left += yy[ik];
            x_left += xx[ik];
        }
        y_left /= ik2;
        x_left /= ik2;
        
        /* down */
        
        y_right = 0;
        x_right = 0;
        for(ik = i + ik1; ik < i + ik1 + ik2; ++ik){
            y_right += yy[ik];
            x_right += xx[ik];
        }
        y_right /= ik2;
        x_right /= ik2;
        
        /* delta */
        delta = y_left - y_right;
        if(delta > max_delta){
            max_delta = delta;
            x0 = (x_right + x_left) * 0.5;
        }
    }
    
    return x0;
}

double SphData::derivate(std::vector<double> & tt, std::vector<double> & xx, std::vector<double> & xt)
{
    int i, n = tt.size();
    
    xt.clear();
    xt.resize(n, 0.0);
    
    for(i = 1; i < n-1; ++i){
        xt[i] = (xx[i+1] - xx[i-1]) / (tt[i + 1] - tt[i-1]);
    }
    
    xt[0] = xt[1];
    xt[n-1] = xt[n-2];
}

double SphData::detectLeft(std::vector<double> & xx, std::vector<double> & yy)
{
    //while();
    
}

double SphData::detectRight(std::vector<double> & xx, std::vector<double> & yy)
{
}

double SphData::detectPeak(std::vector<double> & xx, std::vector<double> & yy)
{
    int n = xx.size();
    double func, max_func;
    double func_x, max_func_x;
    int max_i;
    func = 0;
    
    int i;
    int ik, ik1, sik, w;
    ik1 = 15;
    
    max_func = 0.0;
    max_func_x = 0.0;
    for(i = ik1; i < n - ik1; ++i){
        
        sik = 0;
        func = 0.0;
        func_x = 0.0;
        
        for(ik = i - ik1; ik < i + ik1; ++ik){
            w = ((ik > i)?(i + ik1 - ik):(ik + ik1 - i)) + 2;
            
            sik += w;
            func += yy[ik] * w;
            func_x += xx[ik] * w;
        }
        
        if(sik > 0){
            func = func / sik;
            func_x = func_x / sik;
        }else{
            func = func_x = 0.0;
        }
        
        if(func > max_func){
            max_func = func;
            max_func_x = func_x;
        }
    }
    
    return max_func_x;
}
